function Table3_Figure1

addpath('../Data')
addpath('basic_functions')

Nfactors = 3;
UB = 90;
LB = 10;

[Return,Z,IsF,Nt,Chars,CharsNames,~,~,~] = package_data_all_greeks;

Return_long = Return.ret_long_60_daily;
Return_short = Return.ret_short_60_daily;
Return = Return.ret_daily;

% construct W and X matrices that will be used for IPCA
[N,T,L] = size(Z);
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t;

% run IPCA
[G, F] = IPCA(X,W,Nt,Nfactors);

% compute expected returns with look ahead bias
EReturn = NaN(N,T,'single');
Alpha = NaN(N,T,'single');
Alpha_long = NaN(N,T,'single');
Alpha_short = NaN(N,T,'single');
for t = 2:T
    Betas = squeeze(Z(:,t-1,:))*G;
    EReturn(:,t) = Betas*F(:,t);
    Alpha(:,t) = Return(:,t) - Betas*F(:,t);
    Alpha_long(:,t) = Return_long(:,t) - Betas*F(:,t);
    Alpha_short(:,t) = Return_short(:,t) - Betas*F(:,t);
end

% compute strategies returns
[N,T,L] = size(Chars);
XXX = Chars;
RET = NaN(T,L);
mean_ret = NaN(L,1);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            xx = XXX(~~ind,t-1,i);
            nS = xx<=prctile(unique(xx),90);
            nL = xx>=prctile(unique(xx),10);

            rr = Return(~~ind,t);
            RET(t,i) = mean(rr(nL),'omitnan') - mean(rr(nS),'omitnan');
        end
    end
    mean_ret(i) = mean(RET(:,i),'omitnan');
end

strat_ret = NaN(L,2);
strat_ret_TC = NaN(L,2);
strat_alpha = NaN(L,2);
strat_alpha_TC = NaN(L,2);

% load bootsrap distribution of gammas
load tables_outputs/GammaVariances GammaBoot;
XXX = Chars;
RET = NaN(T,L);
ERET = NaN(T,L);
ALPHA = NaN(T,L);
ALPHA_TC = NaN(T,L);
ZZ = NaN(T,L+1,L);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            rr = Return(~~ind,t);
            err = EReturn(~~ind,t);
            aa = Alpha(~~ind,t);

            rr_long = Return_long(~~ind,t);
            rr_short = Return_short(~~ind,t);

            aa_long = Alpha_long(~~ind,t);
            aa_short = Alpha_short(~~ind,t);
            zz = squeeze(Z(~~ind,t-1,:));
            xx = XXX(~~ind,t-1,i);

            if mean_ret(i) > 0
                nS = xx<=prctile(unique(xx),LB);
                nL = xx>=prctile(unique(xx),UB);
            else
                nS = xx>=prctile(unique(xx),UB);
                nL = xx<=prctile(unique(xx),LB);
            end

            RET(t,i) = mean(rr(nL),'omitnan') - mean(rr(nS),'omitnan');
            RET_TC(t,i) = mean(rr_long(nL),'omitnan') - mean(rr_short(nS),'omitnan');
            ERET(t,i) = mean(err(nL),'omitnan') - mean(err(nS),'omitnan');
            ALPHA(t,i) = mean(aa(nL),'omitnan') - mean(aa(nS),'omitnan');
            ALPHA_TC(t,i) = mean(aa_long(nL),'omitnan') - mean(aa_short(nS),'omitnan');
            ZZ(t,:,i) = mean(zz(nL,:),'omitnan') - mean(zz(nS,:),'omitnan');
        end
    end

    % compute return average and tstat
    ind = isfinite(RET(:,i));
    strat_ret(i,1) = 100*mean(RET(~~ind,i),'omitnan');
    strat_ret(i,2) = sqrt(sum(ind))*mean(RET(~~ind,i),'omitnan')./std(RET(~~ind,i),'omitnan');

    % calculate variance of alpha
    vv = var(ALPHA(~~ind,i));
    for l = 1:L+1
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vv = vv + var(zf)*(var(GammaBoot(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBoot(:,k1,l),GammaBoot(:,k2,l));
                    vv = vv + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_alpha(i,1) = 100*mean(ALPHA(~~ind,i),'omitnan');
    strat_alpha(i,2) = sqrt(sum(ind))*mean(ALPHA(~~ind,i),'omitnan')./sqrt(vv);

    % compute return_TC average and tstat
    strat_ret_TC(i,1) = 100*mean(RET_TC(~~ind,i),'omitnan');
    strat_ret_TC(i,2) = sqrt(sum(ind))*mean(RET_TC(~~ind,i),'omitnan')./std(RET_TC(~~ind,i),'omitnan');

    % compute alpha average and tstat + TC
    vv = var(ALPHA_TC(~~ind,i));
    for l = 1:L
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vv = vv + var(zf)*(var(GammaBoot(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBoot(:,k1,l),GammaBoot(:,k2,l));
                    vv = vv + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_alpha_TC(i,1) = 100*mean(ALPHA_TC(~~ind,i),'omitnan');
    strat_alpha_TC(i,2) = sqrt(sum(ind))*mean(ALPHA_TC(~~ind,i),'omitnan')./sqrt(vv);
end


disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp(' ')
disp(' ')
YY = strat_alpha(:,[1 2]);
bh_t(1,1:2) = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);

X = YY;
for l = 1:size(X,1)/2
    if floor(X(l,2)*1000) >= floor(bh_t(1,2)*1000)
        bb = ['&\textbf{' num2str(X(l,1),'%1.2f') '}&\textbf{(' num2str(X(l,2),'%1.2f') ')}'];
        aa = ['\textbf{' CharsNames{l,1}   '}'];
    else
        bb = ['&' num2str(X(l,1),'%1.2f') '&(' num2str(X(l,2),'%1.2f') ')'];
        aa = CharsNames{l,1};
     end
    if floor(X(size(X,1)/2+l,2)*1000) >= floor(bh_t(1,2)*1000)
        cc = ['&\textbf{' num2str(X(size(X,1)/2+l,1),'%1.2f') '}&\textbf{(' num2str(X(size(X,1)/2+l,2),'%1.2f') ')}'];
        dd = ['\textbf{' CharsNames{size(X,1)/2+l,1} '}'];
    else
        cc = ['&' num2str(X(size(X,1)/2+l,1),'%1.2f') '&(' num2str(X(size(X,1)/2+l,2),'%1.2f') ')'];
        dd = CharsNames{size(X,1)/2+l,1};
    end

    disp([aa bb '&&' dd cc '\\'])
end
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')

figure(1)
set(gcf,'Position',[1 1 750 500])
sig = floor(1000*strat_alpha(:,2)) >= floor(1000*bh_t(1,2));
plot(strat_ret(~sig,1),strat_alpha(~sig,1),'*')
hold on
plot(strat_ret(sig,1),strat_alpha(sig,1),'r*')
set(gca,'ylim',[-1 3])
set(gca,'xlim',[0 3])
xlabel('Raw return')
ylabel('IPCA alpha')
hline = refline(0,0);
hline.LineStyle = ':';
vline = refline(1,0);
vline.LineStyle = ':';
saveas(gcf,'figures_outputs/Figure1.eps','epsc')
saveas(gcf,'figures_outputs/Figure1.jpg','jpg')

return